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Abstract 

With increasingly “big” data available in biomedical research, deriving accurate and reproducible 
biology knowledge from such big data imposes enormous computational challenges. In this paper, mo¬ 
tivated by recently developed stochastic block coordinate algorithms, we propose a highly scalable ran¬ 
domized block coordinate Frank-Wolfe algorithm for convex optimization with general compact convex 
constraints, which has diverse applications in analyzing biomedical data for better understanding cellular 
and disease mechanisms. We focus on implementing the derived stochastic block coordinate algorithm 
to align protein-protein interaction networks for identifying conserved functional pathways based on the 
IsoRank framework. Our derived stochastic block coordinate Frank-Wolfe (SBCFW) algorithm has the 
convergence guarantee and naturally leads to the decreased computational cost (time and space) for 
each iteration. Our experiments for querying conserved functional protein complexes in yeast networks 
confirm the effectiveness of this technique for analyzing large-scale biological networks. 


1 Introduction 

First-order methods in convex optimization have attracted significant attention in statistical learning in re¬ 
cent years. They are appealing to many learning problems, such as LASSO regression and matrix completion, 
which have diverse applications in analyzing large-scale biological systems and high-dimensional biomedical 
measurement profiles m^- These first-order optimization methods scale well with the current “big” data in 
many biomedical applications due to their advantages that they have low computation burden per iteration 
and they are easy to be implemented on parallel computational resources. 

In this paper, we focus on the Frank-Wolfe algorithm, which is also known as the conditional gradient 
method. One of its advantages is that at each iteration step it decomposes the complex constrained opti¬ 
mization problem into sub-problems which are easier to solve. Additionally, it is a projection free algorithm, 
which avoids solving the projection problem for constrained optimization as done in many other algorithms. 
The original Frank-Wolfe algorithm, developed for smooth convex optimization on a polytope, dates back 
to Frank and Wolfe |1]. Dunn and Harshbarger [Hi have generalized the algorithm to solve the optimiza¬ 
tion for more general smooth convex objective functions over bounded convex feasible regions. Recently, 
researchers have proposed stochastic optimization ideas to scale up the original Frank-Wolfe algorithm. 

Based on these previous seminal efforts, our main contribution in this paper is that we generalize the 
stochastic block coordinate Frank-Wolfe algorithm proposed in [^, previously with block separable con¬ 
straints, to solve more general optimization problems with any convex compact constraints, including the 
problems with block inseparable constraints. Such a generalized algorithm has a broader range of biomedical 
applications, including biological network alignment. We prove the convergence of our generalized stochas¬ 
tic block coordinate Frank-Wolfe algorithm and evaluate the algorithm performance for querying conserved 
functional protein complexes in real-world protein-protein interaction (PPI) networks. 

In the following sections, we first describe the model formulation of the optimization problems that 
we are generally interested. Specifically, to address potential difficulty from more general convex compact 
constraints, we derive a new stochastic block coordinate Frank-Wolfe algorithm and provide the convergence 
proof. Then, we formulate the IsoRank problem for network alignment m as a convex programming 
problem and develop a SBCFW-IsoRank algorithm based on our new stochastic block coordinate Frank- 
Wolfe algorithm. At last, in our experiments, we show the efficiency and effectiveness of our algorithm for 
solving the PPI network query problem. 
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2 Stochastic Block Coordinate Descent Prank-Wolfe Algorithm 


Consider the minimization problem; 

min: /(x) 
s.t. X e D, 


( 1 ) 


where the objective function /(x) is convex and differentiable on R^, and the domain 23 is a compact convex 
subset of any vector space. We assume that the optimal solution x* to the above problem is non-empty and 
bounded without loss of generality. 

Assume that we can decompose the solution space into n equal-size subspaces: 


n 

R^ = 0R^-, 7V = y^W, 

7—1 


( 2 ) 


where Ni = ... = Ni = ..., and R'^* denotes the ith equal-size subspace along the corresponding coordi¬ 
nates. This decomposition enables scalable stochastic optimization algorithms. Based on this decomposition, 
we introduce matrices Ui, who sum up to an identity matrix In = X]r=r Ui is a N X N matrix with 

Ui(t, t) = l,t & R^' on its diagonal and the other entries being equal to zero. In typical stochastic optimiza¬ 
tion algorithms |12] . instead of computing the gradient V/(x) at each iteration, the partial gradient of /(x) 
on a randomly selected subspace R^’ is used: 


VJ(x) = C/.V/(x). (3) 

Now we generalize the previous stochastic block coordinate Frank-Wolfe algorithm derived in to solve 
more general optimization problems with any compact convex constraints 23. The new generalized stochastic 
block coordinate Frank-Wolfe (SBCFW) algorithm is illustrated in Algorithm 1. In the pseudo code, the 
operation i = randomly selects one of the n equal-size subspaces to update the partial gradient at each 
iteration with the same probability. In addition, Uj x s = Uj x x*^ denotes the condition that the elements 
of the jth block of s equal to the elements of the jth block of x^. 


Algorithml: Generalized SBCFW Algorithm 

1 Let x*^ s 23, k = 0. 

2 While Stopping criteria not satished do 

3 Randomly divide R^ into n blocks R^ = 0 R^*; 

i=l 

4 Choose i = €k; 

5 Find sf’ such that 

6 := arg min Vi/(x^)^(s — x^); 

Uj xs—Uj xx^ 
sGS 

7 Determine the step size 7 

8 'y := arg min /((I — yjx^ 0 ys^); 

7 e[o.i] 

9 Update x'^+i := (l- 7 )x '=0 7 sf; 

10 k = k + 1; 

11 Endwhile 


Note that our generalized SBCFW algorithm is similar to the algorithm in [9], which aims to solve opti¬ 
mization problems with block separable constraints and has the sub-linear convergence property. However, 
our algorithm provides a more generalized framework, which can manipulate any convex and compact con¬ 
straints no matter whether they are block separable or not. Because the setup of our algorithm is more 
general without any specihc structure, it is difficult to obtain theorectical convergence rate guarantees. In 
this paper, we only provide the proof that our SBCFW converges to the global optimum. The convergence 
guarantee of the generalized SBCFW algorithm is provided by Theorem 1 below, which is based on Lemma 
1 : 
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Lemma 1: At each iteration of the SBCFW algorithm, the following inequality holds 


(Ajsf] - < 0, (4) 

where is the expectation of with respect to the random selection of the Ah cordinate block to the 

corresponding subspace. 

Proof. Assuming at the kt\i iteration, we solve the following optimization problem: 

min: Z^(s) := V,/(x'")'^(s - x'^) 

s.t. Uj X s = Uj X x^, Vj ^ i, (5) 

s e D. 


The solution to (H is s^. With sf achieving the minimum of (H , we have 

< Z^(x'=) = V./(x'=)^(x'= - x'^) = 0. 


( 6 ) 


Therefore, 

ZUs,^)=V./(x'=f(s,^-x'=)<0. 

Taking expectation on both sides of the above inequality with respect to random blocks, we obtain 

A, [V,/(x'=)^(sf - x'^)] < 0 

^ l^V./(x'=f(s,^-x'=) <0 

n 

i 

^ <0 
^ V/(x^'f (A,[s,^]-x'=) <0. 


(7) 


( 8 ) 


The inequality in the third line can be derived based on the fact that ^ — x^ is a vector with only its Ah 
coordinate block having non-zero values and the other parts being all zeros. With that, the summation in 
the second line can be written as the inner product between vectors and ~ ^ 

We now analyze the convergence of the new SBCFW algorithm based on Lemma 1 from two cases. The 
first case is when 

V/lx'^f (F;,[sf]-x'=) =0. (9) 

This simply means that x^ is a stationary point. Because the original objective function /(x) is convex, we 
can conclude that x^ is the global minimum. Another case is when 

v/ix'^f (f;,[s,"]-x'=) <0, (10) 

indicating that Ei[s^] — x*^ is a decent direction based on the definition m- Hence, Ei[s^] — x^ can move 
along the direction to get closer to the global minimum in expectation. Furthermore, we compute the optimal 
step size at each iteration, therefore the objective function values are guaranteed to be non-increasing. With 
that, we present Theorem 1 as follows: 

Theorem 1: The sequence {/(x^),/(x^),...,/(x^), ...} generated by the SBCFW algorithm is non¬ 

increasing 

/(xi) > /(x2) > ... > /(x'^) > /(x'^+i), fc^OO. (11) 
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3 Biological Network Alignment 


3.1 Optimization Model Formulation 

In this section, we re-formulate the involved optimization problem for the network alignment algorithm— 
IsoRank m to address the potential computational challenges of aligning multiple large-scale networks. The 
new formulation has the same mathematical programming structure as the problem Q. 

Let Ga and Gh be two biological networks to align. Two networks has Na and TVf, vertices respectively. 
We define B G Cartesian product network from Ga and Gb- B = Ga®Gb- Denote 

the all one vector 1 G and 

B = B X Diag(Bl)-\ (12) 

where Diag(i31) can be considered as a degree matrix with B1 on its diagonal and all the other entries equal 
to zero. B contains the transition probabilities for the underlying Markov random walk in IsoRank HH. It 
is well known that if Ga and Gb are connected networks and neither of them is bipartite graph, then the 
corresponding Markov chain represented by B is irreducible and ergodic, and there exists a unique stationary 
distribution for the underlying state transition probability matrix B. The goal of the IsoRank algorithm 
is to find a right maximal eigenvector of the matrix B: Bx = x and l^x = 1, x > 0, which corresponds 
to the best correspondence relationships between vertices across two networks. When two networks are 
of reasonable size, spectral methods as well as power methods can be implemented to solve the IsoRank 
problem m- However, with large-scale networks, the transition probability matrix B can be extremely 
large (quadratic with Na x Nb) and spectral methods can be computationally prohibitive. In this paper, we 
re-formulate this problem of searching for maximal right eigenvector as a constrained optimization problem: 

1 11 _ 11 2 
min: fix) := - \\Bx — x 

/ 2 II II (13) 

s.t. l^x = 1, X > 0. (fi) 

After expanding the objective function, we obtain f(x) = -x^Mx, where M = B^B — B — B^ + !■ Therefore 
the equivalent optimization problem is 

min: fix) := -x^Mx 

' 2 (14) 

s.t. l^x = 1, X > 0. (Jo) 


The gradient of /(x) can be easily computed V/(x) = Mx. Furthermore, we find that the Hessian matrix 
of fix) is M, which is a positive semi-definite matrix proven by Lemma 2: 

Lemma 2: M = B^B — B — -|- / is positive semi-definite. 

Proof. M can be written as M = (B — I)'^ (B — I), which proves the lemma. □ 


With Lemma 2, it is obvious that the objective function /(x) is convex. Also, the constraint set 
= {x|x^l = 1, X > 0} is a unit simplex, which is convex and compact. Hence, the IsoRank problem (131 has 
the same problem structure as Q and our generalized SBCFW algorithm can be used to solve (14) with much 
better scalability and efficiency due to the efficiency of the randomized partial gradient computation at each 
iteration. Similarly as in m, in addition to network topology, we can incorporate other information in the 
formulation for more biologically significant alignment results by replacing B with B = aB + il — a)Sl'^, a G 
[0,1]. Here S = S'/] S'] is a normalized similarity vector with size Na x Nb, cancatenated from the doubly 
indexed similarity estimates ^([mju]) based on the sequence or function similarity between vertices u in Ga 
and V in Gb. 


3.2 SBCWF-IsoRank Algorithm 

As shown in Section 3.1 /(x) in © is convex and the constraint set h in ( |14[ ) is a convex compact set. 
Therefore, we can apply the generalized SBCWF algorithm proposed in Sectionj^to solve the corresponding 
optimization problem (14). The detailed algorithm is illustrated in Algorithm 2. Here we want to emphasize 
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that, in each iteration of our SBCFW-IsoRank algorithm, both the time complexity and the space complexity 

f N'^\ 

are O ( — I, which is achieved through tracking the vectors of pk = and = ifsf at step 2 and 10 
\ n J 

of each iteration in Algorithm 2, respectively. The stopping criterion is ||ilx — x|| < ^ ||x||, which can be 
efficiently estimated by 

||i?x — x|| = x^Mx = (£'x)^i?x = p^Pfc, (15) 

which is taken in line 11 in the SBCFW-lsoRank algorithm. 


Algorithm 2: SBCFW-IsoRank Algorithm 

Input: n and E 

1 for k = 0, ..., oo do 

2 randomly divide into n equal-size parts 

3 choose i = €k 

4 if(fc == 0) 

TL 

5 initialize the tth block of x° with — 

N 

6 endif 

7 compute pfc = Ex.'^ and = [E'^]iPk 

8 solve the sub-problem: 

9 s^'-arg min V*/(x'=)^(s - x'=) 

UjXs—UjXx^, 

sGSj, 

10 compute = Esf 

11 if p^pfe < Cljxll 

12 break; 

13 endif 

14 compute the step size 7 ^: 


15 


7fe = 


min{ 7 , 1 } 

0 


7 > 0,7 


p^Pfc -Pfcqfc 

p^Pfc - 2p|’qfc -f ql’qfe 

o.w. 


16 x'^+i =x'=-H 7 '=(s,^-x'=) 

17 endfor 
Output: x'^ 


3.3 Initialization 

In order to guarantee both the time and space complexity to be O(^) at each iteration, we can not initialize 
the algorithm with randomly generated x° to avoid a multiplication of a matrix of size N x N and a vector 
of size N, whose time and space complexity would be 0{N^). We propose to initialize x° in the following 
way: First, randomly divide R'^ into n parts with equal sizes and randomly pick the ith part. Then, 
we initialize every elements in the iih part with which makes x*^ in the feasible space defined by the 

constraint set Using the above initialization strategy, the time and space complexity for computating 
Vi/(x°), po = Ex^ and qo = Es^ are all under O(^), which is easy to verify. 


3.4 Algorithm to Solve the Sub-problem 


As shown in the SBCFW-IsoRank algorithm, at each iteration we need to solve a sub-problem. Fortunately, 
the sub-problem can be solved in a straightforward manner for the optimization problem (14). For the 
following sub-problem at iteration k: 


min: V J{x^f {s - x>^) 

s.t. s G (16) 

Uj X s = Uj X x^, Vj ^ i, 
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the optimal solution is s* = + Lej , where is an all-zero vector except that the jth element is 

1 and L = x^(Z). Here, j is the index of the coordinate with the smallest value in the zth block of 

j = arg min- [Vi (17) 


iGR* 


3.5 Optimal Step Size 

To obtain the optimal step size at each iteration, we need to solve the following optimization problem: 

min: (x'^ + 7 (s'= - x'^))^ M (x'^ + - x'^)) 

s.t. 0 < 7 < 1 , 


(18) 


which is the classic quadratic form with respect to 7 . If 7 = 


p^Pfe -prqfc 


Pfc Pk - 2p^ qfe + qi Qfc 


> 0 , which is the 


solution to (181 without any constraints, the optimal solution 7 * is the minimum value between 1 and 7 , 


otherwise 7 * = 0. The definition of pk and are given in lines 7 and 10 in Algorithm 2. 


3.6 Time and Space Complexity 

At each iteration, the most computationally expensive operations are the updates of and q^ (lines 7 and 
10 of SBCFW-lsoRank) and the calculation of the partial gradient Vi/(x^) (line 7 of SBCFW-lsoRank). 
The calculation of p* and q^ are similar. From line 10 of Algorithm 2, we know 


Pfc = 

= E{x- 

= Pfc-i 


fc-i 


■7 


7' 

fc-i 


,_i(sfc-i_xfc-i)) 

£;(s'=-i-x'=-i). 


(19) 


The second equation is derived by replacing x^ with the equation in line 16 of our SBCFW-lsoRank algorithm. 
Because we keep tracking p^ at each iteration, we do not need to recompute Pfc_i. Therefore, we only need 

/ 

to compute — x^'“^), which takes O - I operations because (s^“^ — x^“^) is a vector, with only 


n 


its Bh block being non-zeros and all the other parts are zeros. Additionally, the memory consumption is also 
O ( - 1 by the similar argument. Similarly, we can compute q^ 


qfc = Es'^ 

= E{k'^ + (Lej - Uix'^)) 

= Pfc + E{Lej - t/jx'"), 


( 20 ) 


where {Lej — Uix'^) is also a vector with only the Ah block having non-zero values. Therefore, the computation 

of qfc also takes O — operations and consumes O — memory. 

\ n J \ n J 

The equation of calculating Vi/(x^) is as follows: 


VJ(x'=) = [E^lPk, 


( 21 ) 


where the operator is to get the rows of the matrix corresponding to the ith coordinate block. Hence, it 
is easy to verify that the time complexity and space complexity of computing Vif {x^) are O — 


In summary, based on the above analyses, both the time complexity and space complexity of our SBCFW- 

' N‘^\ 

IsoRank at each iteration are O ( - . 

n J 
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Figure 1: Query subnetwork and its aligned result in the target network: (a) the query subnetwork in the 
Krogan’s yeast PPI network; (b) the aligned result in Collins yeast PPI network. 


4 Experiments 

In this section, we apply our SBCFW-IsoRank algorithm to two network query problems. For the first set of 
experiments, we take a known protein complex in an archived yeast protein-protein interaciton (PPI) network 
in one database |8] as the query to search for the subnetwork in another yeast PPI network [5] with different 
archived interactions. We call this yeast-yeast network query problem. The goal of this set of experiments 
is to check the correctness of our algorithm as we have the ground truth for the target subnetwork. With 
that, we aim to test the convergence property of our algorithm under different partitions and the relationship 
between the number of iterations and number of partitions. The second experiment is to query a large-scale 
yeast PPI network in IntAct [6] to hnd similar subnetworks of proteins with similar cellular functionalities 
for a known protein complex in human PPI network. The aim of this experiment is to show that our new 
algorithm can help transfer biology knowledge from model organisms to study potential functionalities of 
molecules in different organisms. 


4.1 Yeast-Yeast Network Query Problem 

We test our SBCFW-IsoRank algorithm on the yeast-yeast PPI network query problem by solving the 
optimization problem introduced in the previous section. We take a subnetwork with 6 proteins (Fig. [^a)) 
from the Krogan’s yeast PPI network [5] as the query example to search for the conserved functional complex 
in a target network, which is the Collins network with 1,622 proteins and 9,074 interactions. The query 
subnetwork is the transcription factor TFIIIC complex in the Krogan’s network and we are interested in 
testing whether we can find the same subnetwork in the Collins network. The dimension of our optimization 
problem is 6 x I, 622 = 9, 732. We run this preliminary example so that we can compare our stochastic 
optimization results with the results from the power method, which is typically done in the original IsoRank 
algorithm m- Theoretically, the time and space complexity of our SBCFW-IsoRank at each iteration are 
both 0{N‘^/n) based on the analysis in Section 3.6 Compared to 0{N‘^) time and space complexity for the 


power method by IsoRank HU, our SBCFW-IsoRank algorithm can scale better with properly selected n. 

As both the query example and the target network contain interactions among proteins from the same 
organism—yeast, we can easily check the correctness of the query result. We define the accuracy as the 
number of corrected aligned proteins divided by the total number of proteins in the query subnetwork. We 
implement the SBCFW-IsoRank algorithm for different numbers of partitions n but use the same stopping 
criterion: Bx — x < ^ ||x|j , ^ = 0.1. In Table I, we hnd that our stochastic optimization algorithm obtains 
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Figure 2: The change of the objective function values with the increasing number of iterations with different 
numbers of partitions. 


the same biologically meaningful results as the power method. 

Fig.[3 shows the changes of the objective function values with respect to the increasing number of 
iterations. As illustrated in Fig. our algorithm converges for all different n. Additionally, we find the 
larger the number of partitions n is, the larger the number of iterations we need to have the algorithm 
converge to the global optimum with the same stopping criterion. This clearly demonstrates the tradeoff 
between the efficiency and scalability of the stochastic optimization algorithms. Interestingly, we notice 
that for n = 10, 30, and 50, the number of iterations does not increase much, which indicates that we may 
achieve fast computation with reasonably large n because our algorithm is more efficient for larger n at each 
iteration. 


Table 1: Comparison on different decompositions with ^ = 0.1. 


of partitions (n) 

Computational time (s) 

#. Iterations 

Accuracy 

2 

11.60 

648 

100% 

5 

8.53 

1,045 

100% 

10 

7.44 

1,742 

100% 

30 

5.05 

1,880 

100% 

50 

4.93 

2,070 

100% 

100 

7.06 

2,942 

100% 

200 

13.05 

4,478 

100% 


To further investigate the performance with different n, we run our algorithm 10 times for each n and show 
the average computational time, the average number of iterations, and the average accuracy score in Table 1. 
From Table I, we observe that for all different n, our algorithm can obtain 100% accuracy, which again 
demonstrates the effectiveness and convergence of our generalized SBCFW algorithm. Also, we notice that 
with the increasing n, the number of iterations increase; however, the computational time is first reducing 
then increasing. For example, when n = 2, our algorithm converges with the smallest number of iterations. 


but the computational time is not the best because at each iteration the algorithm takes O 


iV2 

~Y 


operations. 


In contrast, when n = 50, though the number of iterations is larger, but it reaches the global optimum with 
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Figure 3: Querying human protein complex in a yeast PPI network. The proteins are annotated by their 
gene names. The solid lines are protein interactions and the dash lines denote orthologous relationships 
based on protein sequence similarity by BLAST between the proteins in different organisms, (a) Human 
proteasome core complex; (b) The aligned proteasome core complex in yeast found by SBCFW-IsoRank. 


the least computation time, which is indeed twice faster than n = 2. The trend of the computational time 
implies that there may exist the best number of partitions n*. Empirically, when n < n* the computational 
time decreases while the computational time can increase when n > n*. However, it is difficult to provide a 
theoretical proof for this observed phenomenon. Finally, for the scalibility of the algorithm, we always prefer 
larger n to make the memory requirement as low as possible. 

4.2 Human-Yeast Network Query Problem 

We further study the biological signficance of network query results by our SBCFW-IsoRank algorithm. We 
extract a subnetwork as a query example from a human PPI network archived in IntAct [S]. The query 
subnetwork is the proteasome core complex, with induced interactions among the corresponding proteins 
from IntAct [S]. The proteasome core complex in human consists of 14 proteins in total, as shown in Fig.|^a). 
The target network is the yeast PPI network, also obtained from IntAct [B], which has 6,392 proteins and 
77,065 interactions. Our goal is to find the most similar subnetwork to the human proteasome core complex 
in the target yeast PPI network, based on both the interaction topology and the protein sequence similarity, 
which is computed by BLAST |14) . 

We first construct the alignment network, which has TV = 14 x 6, 392 = 89,488 vertices. By our SBCFW- 
IsoRank algorithm with n = 300, instead of operating a matrix of size 89,488 x 89,488 by the power method, 
we only need to handle a matrix of size 298 x 89,488. At each iteration, the computational time as well 
as the memory requirement are reduced 300 times. Our Matlab implementation of SBCFW-IsoRank on 
a MacPro notebook with 8GB RAM takes only around 750 seconds to converge by reaching the stopping 
criteria (^ = 0.1). 

The identified subnetwork in the target yeast PPI network by our algorithm is illustrated in Fig. [^b). 
To evaluate the biological significance of the obtained subnetwork, we check the p-value based on GO (Gene 
Ontology) enrichment analysis using GOTerm Finder [T]. The identified subnetwork is significantly enriched 
in GO term G0:0005839, which is in fact the same proteasome core complex, with p-value 9.552e — 36. This 
experiment demonstrates that our algorithm can find the biologically consistent groups of proteins with the 
same cellular functionalities as the proteins in the query subnetwork, hence with the capability of transferring 
existing biology knowledge in model organisms (yeast for example) to less studied organisms when the group 
of proteins in the query subnetwork require better understanding of their cellular functionalities. 
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5 Conclusions 

In this paper, we generalize the block coordinate Frank-Wolfe algorithm to solve general convex optimization 
problems with any convex and compact constraint set. Our generalized SBCFW algorithm has the conver¬ 
gence guarantee. We re-formulate the IsoRank problem to such a convex programming problem and solve the 
biological network alignment problem by our SBCFW-IsoRank algorithm, which scales better with the size 
of networks under study. The scalability, efficiency, and effectiveness of our algorithm on solving IsoRank are 
demonstrated for real-world PPI network query problems. In our future work, we will consider the derivation 
of the optimal partition number for better tradeoff between computational efficiency and scalability. 
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